library(Seurat)
library(tidyseurat)
library(tidyverse)

demux_sample <- function(path){
  read_delim(path, show_col_types = FALSE) |>
    filter(DROPLET.TYPE == 'SNG') |>
    mutate(.cell = BARCODE, indiv = SNG.BEST.GUESS, .keep = 'none')
}

identify_indiv <- function(data, x){
  data |>
    filter(indiv == x) |>
    count(ind_cov, sort = TRUE) |>
    pull(ind_cov) |>
    pluck(1)
}

join_clues <- function(path, bid){
  freemx <- demux_sample(path)

overlp <- freemx |>
  mutate(.cell = str_remove(.cell, '-.+')) |>
  left_join(clues) |>
  filter(batch_id == bid)

nsample <- overlp$indiv |> unique() |> length()

overlp |>
  ggplot(aes(indiv, group = ind_cov, fill = ind_cov)) + geom_bar() +
  scale_fill_manual(values = DiscretePalette(20))

ind_id <- 0:(nsample-1) |> map_chr(\(x)identify_indiv(overlp,x))

maj <- tibble(indiv = 0:(nsample-1), majority = ind_id) |>
  right_join(overlp) |>
  filter(ind_cov == majority)

message('Fraction of dominant cluster: ',signif(nrow(maj) / nrow(overlp)))
message('Number valid barcodes: ',nrow(maj))
return(maj)}

clues <- read_csv('~/clues1.csv')

batch_read_id <- read_csv(
  'batch_cov,batch_id
  dmx_count_BH7YT2DMXX_YE_0907,1
  dmx_AbFlare-4,2
  dmx_YE_7-13,3')

clues <- clues |>
  select(-batch_id) |>
  left_join(batch_read_id)

clues |> count(batch_id)

clues |> filter(batch_id == 3) |>
  count(ind_cov_batch_cov)

# identify batch id --------
read_tsv('~/append-ssd/nextflowing/scrna-sle-perez2022-p32-34/cellranger/mtx_conversions/CLUES1_POOL03_2/barcodes.tsv', col_names = 'id') |>
  mutate(.cell = str_remove(id, '-.')) |>
  right_join(clues) |>
  filter(!is.na(id)) |>
  count(batch_cov, sort = TRUE)

# do all ind present in all lib in a pool?
dup_bc <- clues |>
  filter(batch_id == 2) |>
  count(.cell, ind_cov) |>
  filter(n > 1)

# do my freemuxlet result match with clues?
## pool1.2 ----------
freemx1.2 <- join_clues('~/append-ssd/nextflowing/scrna-sle-perez2022-3/cellranger/count/CLUES1_POOL01_2/outs/pool1.2.deconv.clust1.samples.gz',1)

## pool1.3 ----------
freemx1.3 <- join_clues('~/append-ssd/nextflowing/scrna-sle-perez2022-3/cellranger/count/CLUES1_POOL01_3/outs/pool1.3.deconv.clust1.samples.gz',1)

freemx1.3 |>
  ggplot(aes(indiv, group = ind_cov, fill = ind_cov)) + geom_bar() +
  scale_fill_manual(values = DiscretePalette(20))

# how freemuxlet umi affect result
freemx1.1 <- join_clues('~/append-ssd/nextflowing/scrna-sle-perez2022-re1/cellranger/count/CLUES1_POOL01_1/outs/pool1.1.deconv.clust1.samples.gz',1)

# how freemuxlet nsample affect result ----------
freemx2.1n19 <- join_clues('~/append-ssd/nextflowing/scrna-sle-perez2022-re1/cellranger/count/CLUES1_POOL02_1/outs/pool2.1n19.deconv.clust1.samples.gz', 2)

## if nsample=16 ?
freemx2.1n16 <- join_clues('~/append-ssd/nextflowing/scrna-sle-perez2022-re1/cellranger/count/CLUES1_POOL02_1/outs/pool2.1.deconv.clust1.samples.gz',2)

# how removing reads without UB or CB affect -------
freemx1.1 <- join_clues('~/append-ssd/nextflowing/scrna-sle-perez2022-re1/cellranger/count/CLUES1_POOL01_1/outs/pool1.1.deconv.clust1.samples.gz', 1)

freemx1.1cb <- join_clues('~/append-ssd/nextflowing/scrna-sle-perez2022-re1/cellranger/count/CLUES1_POOL01_1/outs/pool1.1n16.deconv.clust1.samples.gz', 1)

freemx1.1 |>
  ggplot(aes(indiv, group = ind_cov, fill = ind_cov)) + geom_bar() +
  scale_fill_manual(values = DiscretePalette(20))

g1<- last_plot()

freemx1.1cb |>
  ggplot(aes(indiv, group = ind_cov, fill = ind_cov)) + geom_bar() +
  scale_fill_manual(values = DiscretePalette(20))

g2<- last_plot()

g1 / g2

# cellranger vs starsolo? -------
freemx3.2cr <- join_clues('~/append-ssd/nextflowing/scrna-sle-perez2022-p32-34/cellranger/count/CLUES1_POOL03_2/outs/pool3.2n16.deconv.clust1.samples.gz', 3)

bc3.2ss <- read_tsv('~/append-ssd/nextflowing/scrna-sle-perez2022-sp32-34/star/CLUES1_POOL03_2/CLUES1_POOL03_2.Solo.out/Gene/filtered/barcodes.tsv.gz', col_names = '.cell')

clues |>
  filter(batch_id == 3 & .cell %in% bc3.2ss$.cell) |>
  filter(.cell %in% bc3.2cr) |>
  nrow()

bc3.2cr <- read_tsv('~/append-ssd/nextflowing/scrna-sle-perez2022-p32-34/cellranger/mtx_conversions/CLUES1_POOL03_2/barcodes.tsv', col_names = '.cell')

bc3.2cr <- bc3.2cr$.cell |> str_remove('-.')

clues |>
  filter(batch_id == 3 & .cell %in% bc3.2cr) |>
  nrow()

clues |>
  filter(batch_id == 3 & .cell %in% bc3.2ss$.cell & .cell %in% bc3.2cr) |>
  nrow()

# merge into single rds ------
## pool3
demux_path <-
list.files(path = '~/append-ssd/nextflowing/',
           pattern = 'samples.gz',
           recursive = TRUE,
           full.names = TRUE)

pool123 <- tibble(
  path = demux_path,
  lib_name = demux_path |> str_extract('pool...'),
  pool_id = demux_path |> str_extract('(?<=pool).') |> as.numeric()
)

comp123 <- pool123

meta_comp123 <- map2(comp123$path, comp123$pool_id, join_clues) |>
  map2(comp123$lib_name, \(x,y)mutate(x, lib_name = y)) |>
  list_rbind()

meta_comp123 <- meta_comp123 |>
  mutate(.cell = str_c(lib_name, '_', .cell)) |>
  select(-...4)

meta_comp123 <- meta_comp123 |>
  count(.cell) |>
  filter(n == 1) |>
  left_join(meta_comp123)

## read h5 ---------
dir3_path <-
  list.files(path = '~/append-ssd/nextflowing/',
             pattern = 'filtered_feature_bc_matrix',
             include.dirs = TRUE,
             recursive = TRUE,
             full.names = TRUE) |>
  str_subset('h5', negate = TRUE)

h5.123 <- tibble(
  path = dir3_path,
  lib_name = demux_path |> str_extract('POOL....') |> str_to_lower() |>
    str_replace('_','.') |> str_remove('0'),
  pool_id = demux_path |> str_extract('(?<=POOL0).') |> as.numeric()
)

h5.123 <- h5.123 |>
  filter(str_detect(lib_name, '1.4|2.2'))

dir3_path <- set_names(h5.123$path, h5.123$lib_name)

# 4m46s to read pool1-3
system.time(merged_mtx <- Read10X(dir3_path, strip.suffix = TRUE))

glimpse(merged_mtx)

sobj123 <- CreateSeuratObject(merged_mtx, min.cells = 3, min.features = 200)

sobj_real <- sobj123 |> filter(.cell %in% meta_comp123$.cell)

sobj_real <- sobj_real |> left_join(meta_comp123)

sobj_real <- sobj_real |> select(-c(indiv, majority))

sobj_real |> write_rds('mission/perez-sle-pool1-3.rds')

sobj_real <- read_rds('mission/perez-sle-pool1-3.rds')

SeuratDisk::as.h5Seurat(sobj_real, 'mission/perez-sle-pool1-3.h5seurat')

sobj_real |>
  as.SingleCellExperiment() |>
  zellkonverter::writeH5AD('mission/perez-sle-pool1-3.h5ad')

sobj_real <- SeuratDisk::LoadH5Seurat('mission/perez-sle-pool1-3.h5seurat')
